# This file contains the S-plus code that was used to implement the estimators and functionals in "A Measure of Asymmetry, Statistical papers", (53): 971-985, (2012).
# As a fully working example the code herein can reproduce the simulations presneted on the paper.
# More general use of the estimates is also feasible by supplying the corresponding data set.  Use of the functions is self - explanatory and is further facilitated by the comments at he preamble of each function.

#N.B. Compatibility with R is not guaranteed as the functions were not written for use with R - however, with minimal modifications, their use in the R environment should be feasible.


options(digits=3, echo=FALSE)
remove(objects())
"Epanechnikov"<-
function(x){
	#Epanechikov kernel, square root of 5 is taken as 2.236068 (to save comp. time).
	ifelse(abs(x)< 2.236068, (3/4)*(1-((x^2)/5 ))/2.236068, ifelse(abs(x)>2.236068, 0,0))
}
"IntEpanechnikov"<-function(x)
{
	ifelse(x< -2.236068, 0, ifelse(x> 2.236068, 1, .5- (2.236068*x*(x^2-15))/100))
}

LNormFourthDeriv<-function(x)
{
	dnorm(x) * (6*x^2 - x^4 - 3)
}

"IntKde"<-function(xin, xout, h, kfun)
{
	n<- length(xin)
	arg1<-(sapply(xout, "-", xin))/h	
	arg2<-kfun(arg1)
	arg3<-colSums(arg2)
	arg3/n
}

"Kde"<-function(xin, xout, h, kfun)
{
	n<- length(xin)
	arg1<-(sapply(xout, "-", xin))/h	
	arg2<-kfun(arg1)
	arg3<-colSums(arg2)
	arg3/(n*h)
}



edf<-function(xin, xout)
{
	n<-length(xout)
	nn<-length(xin)
	xin.use<-sort(xin)
	out<-sapply(1:n, function(i, xin.use, xout) length(which(xin.use <= xout[i])), xin.use, xout)
	out/nn
}

DPIbandwidth<-function(xin)
{
	n<-length(xin)
	RofK <- 3* sqrt(5)/25 	
	MuofK <- 3* sqrt(5)	/35
	g<-PilotWidth1(xin)
	arg1<-(sapply(xin, "-", xin))/g
	arg2<- sum(	-LNormFourthDeriv(arg1))/(n^2 * g)
	psiest<- arg2 # use 1 for weib(a=1)
	#cat("\n",RofK / (MuofK^2 * psiest * length(xin)),"\n")
	(RofK / (MuofK^2 * psiest * length(xin)))^0.2
}

etaweib<-function(a)
{
	
	b<-(2*a -1)/a
	b1<-2*(2*a -1)/a

	c<- (3*a-2)/a
	d<- (gamma(c)/(3^c)) - gamma(b) * gamma(b)/(2^b1)
	#d<-ifelse(d<0, -d,d)
	-sqrt(3) * gamma( b) * ( (3^b - 2^(1+b))/(6^b) ) / sqrt(d )#^{-1/2}
}

"SimpsonInt"<- function(xin, h)
{   #xin = data, h=distance of the points at which function is evaluated
    n<-length(xin)
    nn<-1:n
    int0<-xin[1]
    int2n<-xin[n]
    inteven<-xin[seq(3,n-1,by=2)]
    intodd<-xin[seq(2,n-2,by=2)]
    sumodd<-sum(intodd)
    sumeven<-sum(inteven)
    res<-(h/3)*(int0+4*sumodd+2*sumeven+int2n)
    res
}




	
"PilotWidth1"<- function(sample)
{
    #returns Silverman's default width 1.06 A n^(-1/5)
    n<-length(sample)
    StD<-stdev(sample)
    IqR<-diff(quantile(sample, c(.25, .75)))
    tmp<-min(StD, IqR/1.34)
    DefWidth<- 1.06*tmp*n^(-1/5)
    DefWidth
}

RSample<-function(s,dist, p1,p2)
{
	#cat("\n", dist, " ", p1, " ", p2, " ", s,"\n")
	switch(dist, undist = c(runif(s/2)/(0.5+p1)-1, (runif(s/2)-(0.5+p1)+p2*(0.5-p1))/(0.5-p1)), #rweibull(s, shape=p1), 
					weib = rweibull(s, shape=p1),
					lognorm = rlnorm(s),
					norm = rnorm(s),
					uni = runif(s),
					cauchy = rcauchy(s, p1, p2)	,
					fnorm = abs(rnorm(s)	),
					normmixt = c(rnorm(ceiling(p1*s), 0, 1),rnorm(floor((1-p1)*s), 2,2))		
					)

}


FindEtas<-function(s1, s2, dist, p1, p2)
{
sample1<-RSample(s1,dist,p1,p2)
sample2<-RSample(s2,dist,p1,p2)
h1<- PilotWidth1(sample1)#DPIbandwidth(sample1)#
h2<- PilotWidth1(sample2)#DPIbandwidth(sample2) #

Ui1.1<-Kde(sample1, sample1, h1, Epanechnikov)
Ui1.2<-Kde(sample1, -sample1, h1, Epanechnikov)
Ui1<- (Ui1.1+ Ui1.2)/2

Ui2.1 <- Kde(sample2, sample2, h2, Epanechnikov)
Ui2.2 <- Kde(sample2, -sample2, h2, Epanechnikov)
Ui2<-( Ui2.1 + Ui2.2)/2

#Vi1.1<-IntKde(sample1, sample1, h1, IntEpanechnikov)
#Vi1.2<-IntKde(-sample1, -sample1, h1, IntEpanechnikov)
#Vi1<-(Vi1.1+Vi1.2)/2

Vi1 <-IntKde(sample1, sample1, h1, IntEpanechnikov)
Vi2<-IntKde(sample2, sample2, h2, IntEpanechnikov)

Wi1<-edf(sample1, sample1)
Wi2<-edf(sample2, sample2)

MeanUi1<-mean(Ui1 )
MeanVi1<-mean(Vi1)
MeanUi2<-mean(Ui2 )
MeanVi2<-mean(Vi2)
MeanWi1<-mean(Wi1)
MeanWi2<-mean(Wi2)

etahat1<- -(sum(Ui1 * Vi1) - s1 * MeanUi1 * MeanVi1)/sqrt( (sum(Ui1^2) - s1 *MeanUi1^2) * (sum(Vi1^2) - s1 * MeanVi1^2))
etahat2<- -(sum(Ui2 * Vi2) - s2 * MeanUi2 * MeanVi2)/sqrt((sum(Ui2^2) - s2 *MeanUi2^2) * (sum(Vi2^2) - s2 * MeanVi2^2))
etabreve1<- -(sum(Ui1 * Wi1) - s1 * MeanUi1 * MeanWi1)/sqrt((sum(Ui1^2 )- s1 * MeanUi1^2) * (sum(Wi1^2) - s1 * MeanWi1^2))
etabreve2<- -(sum(Ui2 * Wi2) - s2 * mean(Ui2 ) * MeanWi2)/sqrt((sum(Ui2^2) - s2 *mean(Ui2 )^2) * (sum(Wi2^2) - s2 * MeanWi2^2))

etatilde1<- -(sum(Ui1 * Vi1) - s1 * MeanUi1 * 0.5)/sqrt(s1 * (sum(Ui1^2) - s1 *MeanUi1^2) * 1/12 )
etatilde2<- -(sum(Ui2 * Vi2) - s2 * MeanUi2 * 0.5)/sqrt(s2 * (sum(Ui2^2) - s2 *MeanUi2^2) * 1/12 )

c(etahat1, etahat2, etabreve1, etabreve2, etatilde1, etatilde2)
}

Simulation<-function(Iterations,  s1,s2, dist, p1,p2, TrueEta)
{
	matsw1<-matrix(nrow=Iterations, ncol=3)
	matsw2<-matrix(nrow=Iterations, ncol=3)
for(i in 1:Iterations)
{
	alletas<-FindEtas(s1,s2, dist, p1,p2)
	matsw1[i,]<-alletas[c(1,3,5)]#c(etahat1, etabreve1, etatilde1)
	matsw2[i,]<-alletas[c(2,4,6)]#c(etahat2, etabreve2, etatilde2)
}
mses1<-(TrueEta-colMeans(matsw1))^2 +colVars(matsw1)
mses2<-(TrueEta-colMeans(matsw2))^2 +colVars(matsw2)

c(mses1, mses2,  colMeans(matsw1), colMeans(matsw2),  colVars(matsw1), colVars(matsw2))
}
set.seed(20)
final.sim.res<-matrix(ncol=18, nrow = 15)
final.sim.res[1,]<-Simulation(100,30,50, "weib", 2,0, etaweib(2))
final.sim.res[2,]<-Simulation(100,30,50, "weib", 3,0, etaweib(3))
final.sim.res[3,]<-Simulation(100,30,50, "weib", 4,0, etaweib(4))
final.sim.res[4,]<-Simulation(100,30,50, "weib", 5,0, etaweib(5))
final.sim.res[5,]<-Simulation(100,30,50, "fnorm", 1,0, 0.9526)
final.sim.res[6,]<-Simulation(100,30,50, "lognorm", 1,0, 0.909649)
final.sim.res[7,]<-Simulation(100,30,50, "norm", 1,0, 0)
final.sim.res[8,]<-Simulation(100,30,50, "normmixt", 0.945,0, 0.1)
final.sim.res[9,]<-Simulation(100,30,50, "normmixt", 0.101,0, 0.1)
final.sim.res[10,]<-Simulation(100,30,50, "normmixt", 0.872, 0, 0.2)
final.sim.res[11,]<-Simulation(100,30,50, "normmixt", 0.175, 0, 0.2)
final.sim.res[12,]<-Simulation(100,30,50, "normmixt", 0.773,0, 0.3)
final.sim.res[13,]<-Simulation(100,30,50, "normmixt", 0.256,0, 0.3)
final.sim.res[14,]<-Simulation(100,30,50, "normmixt", 0.606, 0, 0.4)
final.sim.res[15,]<-Simulation(100,30,50, "normmixt", 0.382, 0, 0.4)
colnames(final.sim.res)<-c("MseEtaHat30", "MseEtaBr30", "MseEtaTil30", "MseEtaHat50", "MsetaEBr50", "MseEtaTil50", "AvEtaHat30", "AvEtaBr30", "AvEtaTil30", "AvEtaHat50", "AvEtaBr50", "AvEtaTil50", "VEtaHat30", "VEtaBr30", "VEtaTil30", "VEtaHat50", "VEtaBr50", "VEtaTil50")
row.names(final.sim.res)<-c("W 2", "W 3", "W 4", "W 5", "F N", "L N","NOR", "NM1", "NM2", "NM3", "NM4", "NM5", "NM6", "NM7", "NM8")
final.sim.res